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The  requirements  on  SDIO  imaging  problems,  although  somewhat  severe,  can  generally 
be  accomplished  by  using  very  large  optics.  Unfortunately  the  expense  associated  with 
constructing,  testing  and  deployment  is  generally  prohibitive  (as  I  can  attest  having  served  on  the 
committee  investigating  the  Hubble  Space  Telescope  failure)  and  other,  less  expensive 
approaches,  are  desirable. 

The  main  theme  of  this  contract  is  the  development  of  unconventional  imaging 
techniques.  The  various  unconventional  imaging  methods  advocated  result  in  complicated 
nonlinear  relations  between  the  object  and  its  image,  inversion  techniques  are  required  to 
translate  the  measured  data  into  a  meaningful  object.  Although  these  methods  require  a  large 
amount  of  computation,  they  are  relatively  cheap  compared  to  the  costs  of  large  telescopes  which 
rely  upon  physical  size  for  increased  resolution. 

The  technical  approach  is  two-fold:  first,  development  of  the  theory  of  the 
unconventional  imaging  scenario  in  question;  second,  adapt  the  latest  techniques  from  numerical 
analysis  to  carry  out  the  inversion  from  measured  data  to  reconstructed  object.  Any  inversion 
technique  adopted  must  be  robust  with  respect  to  measurement  noise. 

The  contract  was  originally  for  two  years,  but  only  one  year  was  actually  funded.  As  a 
consequence  only  one  of  the  three  research  projects  was  completely  finished.  This  work  has 
been  written  up  for  publication  under  the  title  "Inversion  of  the  Modulus/Phase  of  a  Coherently 
Illuminated  Object  from  its  Measured  Diffraction  Image".  A  copy  of  the  manuscript  is  included 
in  this  section.  Aside  from  its  general  significance,  these  techniques  developed  in  this  paper  are 
of  direct  use  in  optical  microscopy  and  can  be  further  developed  so  as  to  supplement  the  well 
known  technique  of  phase  contrast  microscopy.  The  inversion  algorithm  developed  here  can  also 
be  used  for  lens  testing  and  could,  with  some  further  work,  replace  the  cumbersome  knife  edge 


test  for  wavefront  aberrations. 
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FIGURE  LEGENDS 


Fig.  1  Distribution  of  illuminerice  in  the  dif&action  image  of  a  coherently  illuminated  opaque 
edge;  o(v)  =  0,  u  <  0  and  o(v)  =  1  for  u  >  0  viewed  through  an  annular  aperture 
of  obscuration  radius  €  =  0.05.  The  solid  line  is  for  the  in-focus  situation,  the  dotted 
line  is  for  a  half-wave  of  defocus.  Taken  from  Reference  5. 

Fig.  2  Sample  realization  of  reconstruction  of  modulus  and  phase  of  coherently  illuminated 
object  (dashed  lines)  in  the  presence  of  2%  measurement  “noise”. 

Fig.  3  Sample  realization  of  reconstruction  of  modulus  of  coherently  illuminated  object 
(dashed  line)  in  the  presence  of  4%  measurement  “noise”. 

Fig.  4  Sample  realization  of  reconstruction  of  modulus  of  coherently  illuminated  object 
(dashed  line)  in  the  presence  of  4%  measxirement  “noise”. 
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ABSTRACT 


An  algoritlun  for  recovering  the  modulns  and  phase  of  a  coherently  illununated  object 
from  its  measured  diffraction  image  is  presented.  The  algorithm  is  based  upon  the  fact 
that  both  the  Jacobian  and  Hessian  matrices  can  be  evaluated  exactly  so  that  both  slope 
and  curvature  information  is  available.  The  inversion  problem  is  cast  as  a  nonlinear  uncon¬ 
strained  optimization  problem,  and  trust  region  techniques  are  employed  for  its  solution. 
Representative  numericals  are  presented. 


1.  INTRODUCTION 


The  problem  of  measming  the  diffractiou  image  of  an  incoherently  illuminated  object 
and  working  back  to  determine  (or  estimate)  the  object  itself  has  been  of  great  interest 
in  many  scientific  and  technical  areas  for  the  past  twenty  years  starting  with  the  work  of 
Barakat  and  Blackman  [1]  who  employed  Tichonov  regularization  methods.  The  subject 
of  object  restoration  for  incoherently  illuminated  objects  has  grown  to  such  an  extent  that 
even  listing  the  books  and  major  review  articles  is  a  tedious  undertaking;  however,  we 
should  refer  to  the  standard  work  of  Andrews  and  Hunt  [2]  which  has  served  to  educate 

so  many  people  interested  in  the  topic. 

At  the  other  extreme,  we  have  the  analogous  problem  for  a  coherently  illuminated 
object,  a  much  more  difficult  problem  because  the  relation  between  object  and  diffraction 
image  is  nonlinear,  whereas  the  corresponding  relation  for  incoherently  illuminated  objects 
is  linear.  In  many  respects  the  coherently  iUuminated  situation  is  much  older  than  the 
incoherently  iUuminated  situation  in  that  Ught  microscopists  have  always  encountered 
such  problems.  Their  “solutions”  have  generaUy  been  of  the  old-fashioned  expert  systems 
type;  they  have  encountered  over  the  years  a  variety  of  biological  objects  and  developed 
an  empirical  expertise  in  sorting  out  situations.  In  a  sense  their  primary  artifact  is  the 
diffraction  image,  not  the  actual  object,  as  witness  the  various  phase  contrast  methods 
[3,4].  Unlike  the  incoherent  situation,  there  is  no  guarantee  that  the  topology  of  the  object 
bears  any  resemblance  to  its  coherent  diffiaction  image.  As  an  example  see  Figure  1  which 
shows  the  diffraction  image  of  an  edge  viewed  through  an  optical  system  with  an  annular 

aperture  [5]. 

In  the  present  paper  a  solution  for  the  modulus  and  phase  of  a  coherently  iUuminated 
object  is  obtained  via  a  nonUnear  regularized  minimization.  We  have  brought  to  bear 
powerful  tools  recently  developed  in  numerical  analysis  (particularly  trust  region  consid¬ 
erations)  toward  the  efficient  solution  of  the  inversion  problem.  Even  with  a  CRAY,  the 


Si 


iUustrative  calculations  at  the  end  of  the  paper  required  running  times  of  up  to  several 
minutes. 

The  present  scheme  takes  advantage  of  the  fact  that  we  can  evaluate  both  Jacobian 
and  Hessian  matrices  exactly,  thereby  aUowing  us  to  make  use  of  slope  and  curvature 
information  explicitly.  There  is  no  need  to  make  the  usual  small  residual  approximation 
of  the  Hessian  with  its  convergence  limitation  in  the  presence  of  measurement  noise.  A 
second  benefit  of  knowing  the  Jacobian  and  Hessian  explicitly  is  that  we  can  make  very 
efficient  use  of  the  trust  region  tests  in  determining  the  path  to  local  minima. 

There  are  a  number  of  algorithms  for  the  inversion  of  modulus/phase  of  coherently 
iUuminated  objects  already  published;  see  Stark  [6]  for  an  extremely  useful  summary. 
Practically  all  these  algorithms  will  yield  a  reasonably  accurate  estimate  of  the  support 
of  the  object,  and  the  present  algorithm  is  no  exception.  However,  many  coherently  Ulu- 
minated  objects  contain  changes  in  modular/phase  over  the  surface  so  that  an  algorithm 
(such  as  the  one  disctissed  here)  which  also  can  dehver  information  on  the  distribution  of 
modulus/phase  over  the  object  is  extremely  useful  for  appUcations  (e.g.,  biological  light 

microscopy). 


2.  DIFFRACTION  MODEL 


The  diffraction  image  of  a  coherently  illuminated  object,  I(i,  y)  in  the  image  receiving 
plane  with  coordinates  x,y  is  measured  over  a  square  lattice  of  points  Xm,  Vn- 

xm  =  ^m  rn,n  =  Q,±l,±2,...,±M  (2.1) 

yn  =  Pn 

where  /3  is  a  numerical  constant.  It  is  assumed  that  M  is  large  enough  for 

In  what  follows,  it  will  be  convenient  to  write  (xm,ym)  =  ^mn  as  a  coltimn  vector  I  of 
length  m  =  (2M  +  1)^  using  standard  Fortran  lexiographic  ordering. 

We  assume  that  the  measured  diffraction  image  can  be  modelled  via  scalar  optical 
diffraction  theory  so  that  the  model  diffraction  image,  X(i,y),  is  given  by  the  convolution 

2 

r(x,y)=  J  j  a{x-x',y-y')o{x',y')dx'dy'  (2.3) 

object 

assuming  the  isoplanatic  condition  to  hold. 

The  coherently  illuminated  object  is  characterized  by  a  complex-valued  function 
o{x,y)  with  modulus  lo(a:,y)l  and  phase  arctan[oi(x,y)/o„(z,y)].  The  function  a(x,y) 
is  the  coherent  point-spread  function  of  the  optical  system  performing  the  imaging;  to 
within  multiplicative  factors  it  is 

j  j  A(C,r7)e‘7<*<+3'’'^dC‘iT/  (2.4) 

exit  pupil 


a(i,y)  = 


with,  k  as  the  mean  waveniimber  of  the  coherent  light  an.d  /  the  focal  length.  The  pupil 
function  A(C,  rj)  is  given  by 


^(C,»7)  =  -B(C,?7)e 


ikmc.n) 


with.  W{CiV)  wavefront  aberration  function  over  the  exit  pupil  and  the  ampli¬ 

tude  distribution  over  the  exit  pupil;  both  are  assumed  known  in  the  present  scenario. 

For  the  very  important  case  of  a  circular  aperture  exit  pupil  of  radius  a  for  which 
S(C,  77)  =  1  and  W(C,  77)  =  0,  we  have 


a{x,y)  = 


2/1  [f(i2 +!/")»'"] 


+  y^) 


Since  in  most  applications,  the  exit  pupil  is  circular;  under  these  circumstances  it  is 
convenient  to  employ  normalized  coordinates: 


ka 

u  =  —x, 

C 


V  =  —y 
V 


Consequently  Eqs.  (2.3),  (2.4),  and  (2.6)  become 


X(u,u)  —  J  j  ~  7;')o(u',  v')  du'du' 


lobject 


(2.3a) 


a(u,t,)  =  :  jj  dpJt  f  f 


0<j»^+g*<l 


(2.4a) 


a(u,v)  = 


2Ji((tx2-Hr,2)i/2) 

(«2-1-u2)i/2 


(2.6a) 


Leaving  aside  the  mathematical  details  for  the  moment,  let  us  consider  the  task  before 
us.  We  are  given  measurtd  values  of  the  diffracted  intensity  and  are  reqtured  to  determine 
the  modulus  and  phase  of  the  object  assuming  that  the  measured  diffracted  intensity 
can  be  modeled  as  the  convolution,  Eq.  (2.3),  and  that  we  know  the  parameters  of  the 
optical  system  characterizing  the  coherent  point  spread  function,  Eq.  (2.4).  Said  object 
is  contained  in  the  integrand  of  a  double  integral  which  is  itself  squared.  When  viewed 
from  this  perspective,  we  should  not  be  unduly  optimistic  about  achieving  an  accurate 
solution  because  of  the  inherent  nonlinearity  and  strong  smoothing  action  of  the  double 
integral.  Irrespective  of  the  actual  inversion  method,  the  strong  smoothing  action  of  the 
integration  means  that  much  high  frequency  data  is  lost  and  cannot  be  retrieved,  even  in 
principle.  Only  a  low  frequency  version  of  the  object  can  be  obtained.  Perhaps  the  best 
way  to  consider  the  problem  is  to  interpret  it  thusly:  we  are  given  the  “answer”  (=  effect) 
in  the  form  of  a  noisy  diffraction  image  of  the  object  and  are  attempting  to  determine 
the  “question”  (=  cause),  the  coherently  illuminated  object.  The  mathematical  relation 
between  answer  and  question  is  highly  nonlinear  by  virtue  of  Eq.  (2.2),  bearing  in  mind 
that  a{u,v)  is  additionally  an  oscillating,  complex-valued  function. 

In  a  sense  this  problem  can  be  considered  as  a  two-dimensional  phase  retrieval  problem, 
see  [7-13]  for  various  details.  We  will  not  discuss  the  general  issues  of  ill-posed  inverse 
problems  and  refer  the  reader  to  [14,15]  for  some  of  the  theoretical  aspects. 

Our  diffraction  model  image  I(u,v)  is  to  be  determined  by  the  values  of  the  .image 
at  the  lattice  points  (u^.u^)  inside  a  square  in  the  (u',x;')  plane  that  is  large  enough  to 
contain  the  image.  The  number  of  object  values,  o{v'j,,v't)  =  Okt  is  taken  to  be  N.  It  is 
important  to  remember  that  both  modulus  and  phase  must  be  determined  at  each  of  these 
lattice  points.  Let  n  be  the  number  of  unknown  modulus  and  phase  values  (obviously  n 
must  be  even);  the  •unknown  Okt  axe  to  be  ■written  as  a  column  vector  o. 


In  this  notation  Eq.  (2.3a)  becomes  I(o)  where  J  is  a  vector  of  length  m.  The 
inversion  problem  relates  the  imknown  (complex-valued)  object  values  o  of  the  assumed 
diffraction  model  to  the  measured  diffraction  image  I  via  the  nonlinear  relation 


I=J(o)  .  '  (2.7) 

This  equation  is  to  be  interpreted  as  a  system  of  m  nonlinear  equations  in  h  unknowns. 
Enough  data  must  be  given  to  allow  some  smoothing  of  the  measxired  diffraction  image 
data  as  regards  the  diffraction  model;  consequently  we  let  m  >  n  so  that  the  nonlinear 
system  is  overdetermined.  The  problem  now  reduces  to  the  solution  of  Eq.  (2.7)  in  some 

normed  sense. 

It  is  of  some  interest  to  contrast  the  differences  between  the  incoherently  illuminated 
and  coherently  illuminated  object  situations.  The  relation  between  object  and  image  for 
incoherently  illuminated  objects,  again  assuming  that  the  isoplanatic  conditions  holds,  is 

Ii{x,y)  =  j  J  t{x-x\y-y')oi{x',y')dx'dy'  (2.8) 

object 


where 


Oi(x,  y)  =  intensity  distribution  over  the 
incoherently  radiating  object 


f(x,y)  a  la(x,y)|^  =  incoherent  point-spread  function 


Ii{x,y)  =  intensity  distribution  over  the  difEraction 
image  of  Oi{x,  y) 


Two  points  to  note:  (1)  all  functions  in  Eq.  (2.8)  are  real  and  nonnegative,  (2)  the  relation 
between  answer,  /i(x,  y),  and  question,  Oi(x,  y),  is  Hnear.  On  the  other  hand,  for  the  coher¬ 
ently  iUuminated  object  only  the  diffracted  intensity  is  real  and  nonnegative.  Furthermore 
the  relation  between  I(x,y)  and  o(x,y)  is  nonHnear.  Consequently  the  incoherently  iUu¬ 
minated  object  scenario  is  much  easier  to  invert  than  the  coherently  iUuminated  object 


scenario. 


3.  PRELIMINARIES 


In  order  to  proceed,  we  next  discretize  the  double  integral  on  the  right-hand  side  of 
Eq.  (2.3a) 

.2 


^m)  —  ^mn 


N  N 


m  Ufc,u„  -  vt)o{vk,vt) 


fc=l  ^=1 


(3.1) 


Here  Uk,vt  are  the  quadrature  points  and  ak,at  are  the  corresponding  weight  factors. 

In  the  numerical  algorithm  we  employ  for  the  inversion  (see  next  section),  we  require 
the  first  and  second  derivatives  of  Imn  with  respect  to  Okt  in  order  to  form  Jacobian 
and  Hessian  matrices.  We  need  not  write  out  the  explicit  formulas  because  we  employed 
a  symbol  manipulation  program  to  evaluate  them  directly  in  the  computer  from  which 
numerical  values  are  obtained  internally. 


4.  OUTLINE  OF  SOLUTION  APPROACH 

We  now  outline  the  general  features  of  our  numerical  approach  to  the  nonlinear  phase 
retrieval  problem  via  regularized  unconstrained  minimization.  Our  version  of  the  uncon¬ 
strained  minimization  problem  is  that  of  finding  the  least  value  of  an  objective  function, 
e(o).  The  term  unconstrained  indicates  that  the  variables  o  are  not  limited  in  any  way.  In 
our  application,  we  wish  to  determine  a  global  minimum  of  e,  i.e.,  a  point  o*  satisfying 

e(o)>e(o*),  Vo  .  (4.1) 

Unfortunately  we  must  be  content  with  a  local  minimum: 

e(o)  >  e(o*)  Vo  ,  in  a  neighborhood  of  o*  .  (4.2) 

Solution  of  the  global  minimum  problem  is  far  harder  than  the  solution  of  the  local  mini¬ 
mum  problem  which  itself  is  extremely  complicated  [16,17]. 

The  derivation  of  nearly  all  methods  for  tmconstrained  minimization  is  founded  on  the 
assumption  that  the  objective  function  €(o)  can  be  approximated  by  a  quadratic  function 
in  the  neighborhood  of  a  minimum  point.  Thus  methods  are  sought  which  efficiently 
minimize  quadratic  functions  in  the  hope  that  they  wiU  also  be  effective  on  more  general 
functions,  at  least  in  the  neighborhood  of  a  minimum.  When  first  and  second  derivatives  of 
e(o)  are  known,  such  as  in  our  case,  then  we  can  make  use  of  both  gradient  and  curvatme 
information  to  effect  a  solution.  Furthermore,  both  gradient  and  curvature,  as  governed 
by  G  and  H,  see  Eqs.  (4.4)  and  (4.5),  can  be  calculated  exactly  in  the  context  of  the 
discretized  version  of  our  problem.  Thus,  we  can  avoid  many  of  the  difficulties  associated 
with  situations  where  G  and  H  axe  known  only  approximately.  The  Taylor  series  expansion 
of  the  objective  function  can  be  used  to  approximate  its  minimum  value  from  points  o  near 
to  the  minimum  o*  by  moving  in  a  direction  Ao 


e(o  +  Ao)  «  e(o)  -F  G^(o)(Ao)  -1-  i(Ao)^H(o)(Ao) 


(4.3) 


G  and  H  axe  the  gradient  vector  and  Hessian  matrix  of  the  objective  function,  respectively 
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Note  that  H  is  symmetric.  Both  G  and  H  are  exact  within  the  context  of  the  discretization. 
The  strategy  is  to  determine  the  vector  Ao  of  the  movements  required  to  approximate  the 
rninimnm  from  the  current  point  o. 


Before  proceeding  further,  it  is  important  to  state  the  conditions  for  a  given  point  o  to 
be  an  unconstrained  strong  minimum  o*  (i.e.,  a  point  o*  for  which  the  objective  function 
increases  locally  in  all  directions.  The  first  order  necessary  condition  is  [16,17] 


G(o*)  =  0  .  (4-6) 

However,  this  condition  is  not  sufficient  as  other  types  of  minima  (stationary  points)  also 
satisfy  this  condition.  The  second  order  condition  follows  from  Eq.  (4.2)  and  is  [11,12] 

(Ao)^H(o*)(Ao)  >  0  (4.7) 

which  is  sufficient  to  ensure  that  e(o*  +  Ao)  >  €(o*).  If  e(o)  ^  0,  this  impHes  that  H(o*) 
is  a  positive  definite  matrix.  Therefore  the  second-order  sufficient  condition  for  a  strong 
mini-miiTn  is  that  H{o*)  should  be  positive  defimte. 

Returning  to  Eq.  (4.3),  hereafter  termed  the  local  model  of  the  object  function,  we 
initiate  a  search  for  o*  by  moving  Ao.  This  involves  an  iterative  procedure  in  that  we 


start  from  some  point  o  and  choose  in  some  fashion  a  direction  Ao  in  which  the  minimum 
is  assumed  to  lie.  This  is  repeated  until  the  minimum  is  achieved  (if  possible).  Because  we 
are  employing  a  quadratic  version  of  €(o)  then  the  local  model  of  e(o)  always  allows  us  to 
find  a  solution.  However,  the  local  model  of  €(o)  is  certainly  not  a  useful  approximation 
to  e(o)  itself  except  near  a  minimum.  It  is  an  act  of  faith  that  the  local  model  and  the 
actual  model  are  “close”  in  the  vicinity  of  the  minimum.  Simple  calculations  [16,17]  show 
that  at  a  minimum 

Ao  =  -[H(o)]-1G(o)  .  (4.8) 

This  equation  represents  the  appropriate  direction  Ao  to  take  to  the  minimum  o*,  based 
upon  information  at  o.  This  equation  is  fundamental  to  all  second-order  mimmization 
algorithms  (i.e.,  algorithms  employing  both  G  and  H).  However,  on  the  actual  objective 
function  stirface,  such  as  we  are  faced  with,  the  local  model  of  g(o)  is  only  accurate  in 
the  immediate  vicinity  of  the  minimum.  This  means  that  H  may  not  be  positive  defimte, 
as  required  by  Eq.  (4.7),  since  it  is  evaluated  at  points  other  than  the  minimum.  This 
situation  is  most  Ukely  to  occur  at  some  distance  from  the  minimum  (i.e.,  for  initial  iter¬ 
ations)  since  at  a  point  close  to  the  minimum  all  sufficiently  differentiable  functions  tend 
to  behave  as  a  quadratic  function  as  their  third  and  higher  order  derivatives  in  the  Taylor 
series  become  negUgible.  We  wiU  return  to  the  necessity  of  keeping  H  positive  definite 

very  shortly. 

The  classic  approach  to  limiting  step  size  during  the  iteration  is  a  line  search  [16,17]. 
In  line  search  tactics  we  compute  a  descent  direction  and  subject  this  direction  (which  is 
generally  not  toward  the  unconstrained  minimum)  to  a  minimization  procedure  of  which 
details  can  be  found  in  the  above  references.  Should  the  descent  direction  satisfy  these 
criteria,  this  iteration  is  then  terminated.  We  then  repeat  the  process,  etc.  The  difficulty 
of  practical  implementation  is  two-fold.  First,  such  calculations  are  prohibitively  expen¬ 
sive  and  the  minimization  criteria  are  therefore  only  approximate  to  save  computer  time; 


yet  they  must  be  made  precise  enough  to  ensure  reasonably  quick  convergence.  Second, 
fulfilling  these  contradictory  goals  is  something  of  a  black  art  and  the  programming  effort 
is  very  substantial,  often  occupying  up  to  two-thirds  of  the  coding  for  the  entire  optimiza¬ 
tion.  For  small  problems,  such  as  the  sHt  aperture,  Une  search  methods  are  very  useful 
and  were  employed  along  with  regularized  singular  value  decomposition  in  [18]. 

Although  it  is  possible  to  use  line  search  methodology  for  the  2-D  aperture,  we  have 
chosen  to  employ  a  relatively  new  method  which  offers  many  advantages  over  the  Une  search 
methodology,  the  trust  region  tactic  [16.7],  but  see  also  [18].  A  particular  advantage  of  the 
trust  region  approach  is  in  the  very  strong  global  convergence  properties  which  hold  with 
no  significant  restrictions  on  the  class  of  problems  to  which  they  apply  making  it  especially 
useful  for  the  nonlinear  minimization  of  phase  retrieval  in  two  dimensions.  In  addition,  the 
trust  region  philosophy  requires  less  computation  than  does  the  line  search  philosophy  in 
terms  of  gradient  and  Hessian,  but  more  computations  have  to  be  performed  on  the  local 
model  of  the  objective  function. 

Suppose  we  are  at  point  after  the  fe-th  iteration.  Now  let  us  assume  that  there 
is  a  region  which  we  take  to  be  in  the  shape  of  a  sphere  of  radius 

=  (4.9) 

in  wHch  the  local  model  of  6(o),  given  by  the  Taylor  series  Eq.  (4.2)  a^ees  with  the  actual 
objective  function  in  some  sense.  One’s  obvious  choice  is  to  let 

o(fc+i)  =  o<*’>  -b  .  (4.10) 

where  the  correction  A^**)  minimizes  the  local  model  (A)  for  all  values  of  -}-  A  in  A<*>. 
An  iteration  step  in  o.is  to  be  restricted  by  the  region  of  vaUdity  of  the  Taylor  series.  Thus 
we  compute  the  gradient  and  Hessian,  G  and  H,  appropriate  to  and  minimize  the 


local  model  of  e(o)  in  order  to  determine  the  radius  of  the  trust  region.  In  formal 
language,  we  seek  the  solution  of  the  problem 

minimize  (local  model)  subject  to  constraint  1IA||2  <  .  (4.11) 


Before  we  can  solve  this  subp'roblem  it  is  necessary  to  choose  some  reasonable  criterion 
on  the  radius  of  the  trust  region.  Obviously,  the  criterion  should  not  present  undue 
restrictions,  so  that  should  be  as  large  as  possible  subject  to  some  agreed-upon  idea 
as  to  the  degree  of  closeness  of  the  local  model  of  e(o)  and  e(o).  One  way  is  to  define  the 

quality  coefElcient 


_  .(oW+AW)-6(o(>=))  (4  12) 

6,(oW-bAW)-e,(o(^)) 

where  the  denominator  is  the  predicted  reduction  of  the  local  model  of  e,  now  call  it 
while  the  numerator  represents  the  reduction  in  the  actual  objective  function.  The  closer 
Tjfe  is  to  unity,  the  better  the  agreement  between  e  and  e,.  As  a  stopping  criterion  stop  if 
Tk  >  .8  and  go  to  the  next  iteration.  K  rjk  <  .8  reduce  and  repeat  the  calculation. 
The  Hterature  contains  other  stopping  criteria  but  the  one  quoted  above  seems  adequate 
for  our  purposes.  See  the  quoted  references  for  more  details. 


Thus  far  we  have  outUned  the  general  features  of  the  inversion  problem,  estimating 
o  via  minimization  of  an  objective  function  e,  as  yet  unspecified.  We  feel  that  there  are 
two  objective  functions  of  possible  interest.  The  first  is  the  usual  least  squares  (L2  norm) 

objective  function 

<o)  =  -  lif  .  (4.13) 

t=l 


In  this  version  of  the  problem,  we  consider  the  estimation  of  o  via  an  (unconstrained) 
nonlinear  least  squares  minimization.  A  second  objective  function  is  for  an  Ii  norm 

e(o)  =  ^2 
»=i 


(4.14) 


Consequently  we  have  allowed  considerable  flexibility  in  the  e  minimization  algorithm  so 
as  to  accommodate  both  objective  functions  as  well  as  others  that  may  arise.  For  the 
ptorposes  of  the  present  paper,  we  confine  the  discussion  to  the  I2  norm  aspects. 


Upon  defining  the  vector 


^(o)=  T-I 


which  is  of  length  m,  we  can  rewrite  Eq.  (4.13)  as 


e(o)  =  $^(0)^(0)  . 


(4.15) 


(4.16) 


The  elements  of  G  and  H  can  be  expressed  directly  in  terms  of  the  elements  $  by 
introducing  an  ancillary  matrix  J  given  by 


J  = 


UoT 


This  matrix  is  generally  rectangular.  Thus 


or 


doj 

J  >=i  ^ 


G  =  . 


(4.17) 


(4.18) 


(4.19) 


The  corresponding  elements  of  H  are 

4.0^  6- 

dokdot  ~  ^  dok  dot  ^  *  dokdot 
1=1  1=1 


(4.20) 
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with  =  1,2, . . .  ,n.  The  first  term  on  the  right-hand  side  is  2J^J  and  we  will  call  the 
second  term  S;  thus 

H  =  2J^J  -h  S  (4.21) 

There  is  no  guarantee  that  H  will  always  be  positive  definite  through  the  calculations; 
in  fact  H  will  become  almost  singular  as  we  approach  the  strong  minimum.  To  avoid 
this  state  of  affairs,  we  consider  a  regularized  version  of  H.  In  the  regularized  version, 
we  replace  H  by  H  4-  ?!  where  q,  the  regularizing  parameter,  is  a  small  nonnegative  real 
number.  This  procedure  can  remove  H  from  near  singularity  and  for  sufficient  large  q 
restore  the  positive  definite  character  of  the  Hessian.  In  linear  problems,  there  is  a  well 
established  method  to  estimate  q  [14].  In  the  nonlinear  problem  we  employed  q  in  the 
range  0.01  <  q  <  0.05.  There  was  not  much  difference  in  the  final  answers  as  long  as 
5  >  0;  but  setting  q  =  Q  caused  considerable  numerical  instability  as  expected. 


5.  A  NUMERICAL  EXAMPLE 


It  is  not  our  intention  to  present  a  catalog  of  numerical  results  at  this  time,  yet  we  wish 
to  discuss  the  numerical  workings  of  the  algorithm  in  the  presence  of  “measurement”  noise. 
To  this  end  we  considered  a  known  object  o(u,v)  and  from  it  calculated  the  diffraction 
image  I(u,v)  from  Eq.  (2.3a)  using  the  point-spread  function  a(u,t;)  corresponding  to  an 
in-focus,  aberration-free  optical  system  as  given  by  Eq.  (2.6a).  The  measured  diffraction 
image  I(u,v)  was  taken  to  be  given  by 

I(u,  u)  =  [1  -f-  S^(u,  u)]I(u,  v)  (5-1) 

where  5  is  a  positive  constant  less  than  unity  and  Ur»)  is  a  random  variable  uniformly 

distributed  over  (— 1,  -bl): 

/(m)  =.^  ,  1/^1  <  1 

=  0  ,  (5-2) 

Values  of  6  used  in  the  present  calculations,  such  as  5  =  0.04  are  described  loosely  as 
4%  a  noise.  Note  that  the  noise  we  axe  introducing  is  intensity  dependent  noise;  it  is  not 
intended  to  faithfully  simulate  actual  detector  noise  but  rather  to  mimic  such  noise  for  the 
purpose  of  testing  the  robustness  of  the  inversion  algorithm. 

We  have  chosen  the  following  lone  object  to  illustrate  the  calculations: 

o(u,t;)  =  0,  > 

=  |[3 -h  erf(v)] , 

=  0, 

where  erf(T;)  is  the  error  function.  There  are  two  reasons  for  choosing  this  object.  The 
first  reason  is  that  the  modulus  of  the  object  exhibits  a  spatial  variation.  One  of  the 


—  oo  <  u  <  —15 

— 15  <  u  <  15  (5.3) 

— 15  <  u  <  oo 


(as  yet  unstated)  reqmrements  on  the  algorithm  is  that  it  be  able  to  recover  reasonably 
accTirate  estimates  of  the  object  photometry  in  addition  to  estimating  its  size;  as  noted  in 
the  introduction  almost  any  of  the  currently  available  algorithms  can  perform  this;  very 
few  seem  to  be  able  to  provide  reasonably  accurate  estimates  of  photometry.  The  second 
reason  for  choosing  this  object  is  that  it  is  very  small  in  width  being  slightly  less  than  Airy 
disk  radius  across  (recall  that  an  Airy  disk  radius  «  3.82).  These  two  constraints  present 
a  severe  test  of  the  algorithm. 

Calculations  were  run  for  the  completely  unrealistic  case  of  the  noise-free,  measured 
diffraction  image  of  the  object  given  by  Eq.  (5.3).  The  inversion  algorithm  was. able  to 
retinm  reasonable  answers  as  compared  to  the  true  values.  We  will  not  quote  any  of  these 
results  and  go  to  the  “noisy”  meeisurement  situation. 

In  Figure  2  we  show  a  sample  realization  of  the  reconstruction  in  the  presence  of  2% 
measurement  noise.  Note  the  presence  of  negative  values  of  the  modulus  at  the  edges  of 
the  object;  however  the  photometry  of  the  reconstruction  follows  the  true  object  very  well 
except  at  the  edges.  Two  sample  realizations  of  the  reconstruction  of  the  modulus  of  the 
object,  in  the  presence  of  4%  measurement  noise  are  shown  in  Figs.  3  and  4.  Again  note 
the  unphysical  negative  values  of  the  modulus  at  the  edges  of  the  object.  Considering 
the  fact  that  the  object  is  so  small,  the  modulus  is  in  reasonable  agreement  with  the  true 
object  in  spite  of  measurement  noise. 

The  troublesome  feature  of  the  present  inversion  algorithm  is  obviously  the  occurrence 
of  negative  modulus  values  at  the  edges  of  the  object.  The  edges  axe  somewhat  unrealistic 
because  they  are  of  infinite  slope  and  the  algorithm  tends  to  overshoot  in  the  manner  of  a 
Gibbs  phenomenon.  We  could,  of  course,  consider  objects  with  finite  slopes  at  the  edges 
to  minimize  the  overshoots  and  vmdershoots;  instead  we  have  decided  to  develop  a  vari¬ 
ant  of  the  inversion  algorithm  using  constrained  minimization  to  surmount  this  axmoying 
problem.  Details  will  be  reported  in  the  near  future. 


Although  we  have  run  several  dozen  inversions,  we  do  not  possess  a  sufficiently  large 

base  to  estimate  benchmark  performance.  However,  many  of  our  inversions,  regardless 

/ 

of  the  initial  guess,  took  150-180  iterates  to  converge.  In  a  few  cases,  convergence  was 
achieved  in  half  this  number;  nevertheless  we  also  encountered  some  cases  where  conver¬ 
gence  required  200  iterations. 

Generally  the  algorithm  performance  decomposed  into  two  stages:  (a)  the  initial  stage 
wherein  the  objective  functions  went  through  wild  gyrations  in  magnitude  as  the  trust  re¬ 
gion  subalgorithm  had  to  compute  many  new  gradient  and  Hessian  matrices  while  searching 
for  an  appropriate  direction  of  descent;  (b)  the  second  stage,  when  the  trust  region  subal¬ 
gorithm  has  locked  into  a  subset  of  “optimal”  directions  of  descent,  the  objective  function 
then  undergoes  a  monotone  decrease  to  zero  as  the  number  of  iterations  is  increased.  We 
again  warn  the  reader  that  although  convergence  is  attained,  there  is  no  guarantee  that 
the  global  minimum  has  been  achieved. 

As  a  final  remark,  we  note  that  as  the  analysis  is  heavily  dependent  upon  linear  algebra 
(matrices,  vectors);  it  would  be  of  some  interest  to  consider  doing  calculations  via  parallel 
processing  algorithms  to  speed  up  the  computations  even  more. 


6.  SUMMARY 


Given  the  rather  discursive  presentation,  we  feel  it  is  useful  to  summarize  the  essentials 
of  our  approach  to  the  problem.  We  consider  the  problem  of  determining  the  modulus  and 
phase  of  the  object  at  the  lattice  points  as  one  of  unconstrained  minimization.  A  local 
(i.e.,  quadratic)  model  approach  is  used.  We  are  in  the  fortunate  position  of  employing 
both  the  Jacobian  and  Hessian,  which,  in  the  context  of  the  discretized  diffraction  model, 
can  be  evaluated  analytically!  Thus  we  mahe  use  of  both  slope  and  curvature  information, 
other  methods  only  \ise  slope  information.  The  relatively  new,  and  powerful  trust  region 
algorithm  is  then  utilized  in  place  of  the  usual  line  search  algorithm.  The  trust  region 
algorithm  makes  full  use  of  the  local  model  and  takes  into  account  local  validity;  moreover 
the  algorithm  exhibits  global  convergence  properties. 
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Fig.  1  Distribution  of  illuminence  in  the  diffraction  image  of  a  coherently  illuminated  opaque 
edge:  o(v)  =  0,  v  <  Q  and  o{v)  =  1  for  u  >  0  viewed  through  an  annular  aperture 
of  obscuration  radius  e  =  0.05.  The  solid  line  is  for  the  in-focus  situation,  the  dotted 
line  is  for  a  half-wave  of  defocus.  Taken  from  Reference  5. 


F.g.  2  Sample  realizalioa  of  reconatructioa  of  modulue  aad  phaae  of  cohereatly  iUummated 
object  (daahed  linee)  in  the  presence  of  2%  measurement  “noise”. 
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Fig.  4  Sample  realization  of  reconstruction  of  modulus  of  coherently  illuminated  object 
(dashed  line)  in  the  presence  of  4%  me2isurement  “noise” . 


3  3. 


REPORT  OF  INVENTIONS 


There  were  no  inventions  developed  under  the  aegi 


of  this  contract. 


